#!/usr/bin/tcl


set N [expr {int(pow(2,10))}]
set s_out "synthesized_signal.csv"
set p_out "polar.csv"
set do_x 1





set kmax [expr {int($N / 2)} ]

set pi 3.1415926536

proc fx x {
  global N
  if { $x > [expr { ( $N * 5.0 ) / 16.0 - 1 } ] && $x <= [expr { $N / 16 * 9 - 1} ] } {
    return 1
  } else {
    return 0
  }
}


puts "Calculating basic function."
set s ""
set c ""
for { set k 0 } { $k <= $kmax } {incr k} {
  set ck ""
  set sk ""
  for { set i 0 } { $i < $N } {incr i } {
    set cki [expr {cos( 2 * $pi * $k * $i / $N )}]
    set ski [expr {sin( 2 * $pi * $k * $i / $N )}]
    lappend ck $cki
    lappend sk $ski
  }
  lappend c $ck
  lappend s $sk
}

puts "Calculating ReX and ImX."
set ReX ""
set ImX ""
for { set k 0 } { $k <= $kmax } {incr k} {
  set ReXk 0
  set ImXk 0
  for { set i 0 } { $i < $N } {incr i } {
    set ReXk [expr { $ReXk + [fx $i] * [lindex [lindex $c $k] $i] } ]
    set ImXk [expr { $ImXk + [fx $i] * [lindex [lindex $s $k] $i] } ]
  }
  lappend ReX $ReXk
  lappend ImX $ImXk
}

if { $do_x == 1 } {

puts "Normalizing ReX_ and ImX_."
set ReX_ ""
set ImX_ ""
for { set k 0 } { $k <= $kmax } {incr k} {
  set ReX_k [expr {[ lindex $ReX $k ] / ( $N / 2 )} ]
  set ImX_k [expr {[ lindex $ImX $k ] / ( $N / 2 )} ]
  if { $k == 0 || $k == [ expr { $N / 2 } ] } {
    set ReX_k [ expr { $ReX_k / 2 } ]
  }
  lappend ReX_ $ReX_k
  lappend ImX_ $ImX_k
}

puts "Calculating x_."
set x_ ""
for { set i 0 } { $i <= [expr {$N-1}] } {set i [expr { $i + 1.0 / 32.0} ] } {
  set xi_ 0
  for { set k 0 } { $k <= $kmax } {incr k} {
    set o_rex [lindex $ReX_ $k]
    set o_imx [lindex $ImX_ $k]
    set o_cos [expr {cos( 2 * $pi * $k * $i / $N )}]
    set o_sin [expr {sin( 2 * $pi * $k * $i / $N )}]
    set xi_ [ expr { $xi_ + $o_rex * $o_cos + $o_imx * $o_sin } ]
  }
  lappend x_ $xi_
}


set OFH [open $s_out "w"]
set i 0
foreach lx  $x_ {
  puts $OFH [format "%.6f\t%.6f" $i $lx]
  set i [ expr { $i + 1.0 / 32.0 } ]
}
close $OFH

}

puts "Calculating MagX."

for { set k 0 } { $k <= $kmax } {incr k} {
  set rex [lindex $ReX $k]
  set imx [lindex $ImX $k]
  lappend MagX [expr {sqrt(pow($rex,2)+pow($imx,2))}]
}

set OFH [open $p_out "w"]
set i 0
foreach z $MagX {
  puts $OFH [format "%.6f\t%.6f" $i $z]
  incr i
}
close $OFH

#for { set k 0 } { $k <= $kmax } {incr k} {
#  puts $k
#}

